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ABSTRACT 


This  paper  considers  the  problem  of  estimating  a  two-dimensional 
isotropic  random  field  given  some  noisy  observations  of  this  field  over 
a  disk  of  finite  radius.  By  expanding  the  field  and  observations  in 
Fourier  series,  and  exploiting  the  covariance  structure  of  the  resulting 
Fourier  coefficient  processes,  some  vibrating  membrane  equations  are  obtained 
for  estimating  the  random  field.  These  equations  provide  a  set  of  recursions 
for  constructing  the  field  estimates  as  the  radius  of  the  observation  disk 
increases.  In  the  spectral  domain,  these  recursions  take  the  form  of  Schrodinger 
equations  which  can  be  viewed  as  being  associated  to  an  inverse  scattering 
problem.  , 


Introduction 


In  [1]  -  [5]  seme  vibrating  string  equations  were  obtained  to  estimate 
a  stationary  stochastic  process  given  some  observations  of  this  process  over 
a  finite  interval.  It  was  also  shown  in  [6]  -  [8]  (see  also  [5])  that  the 
estimation  problem  over  a  finite  interval  can  be  formulated  as  an  inverse 
scattering  problem  similar  to  those  appearing  in  the  study  of  transmission 
lines,  or  in  quantum  mechanics.  In  this  paper,  we  will  consider  the  problem 
of  estimating  a  two-dimensional  isotropic  random  field  given  some  observations 
of  this  field  over  a  disk  of  finite  radius.  By  expanding  the  field  and 
observations  in  Fourier  series  and  by  exploiting  the  covariance  structure  of 
the  resulting  coefficient  processes,  some  vibrating  membrane  equations  will  be 
derived  for  solving  a  set  of  filtering  problems  associated  to  the  coefficient 
processes.  These  equations  will  be  used  to  estimate  the  random  field  and  to 
obtain  some  recursions  for  the  optimum  filter  as  the  radius  of  the  obser¬ 
vation  disk  increases. 

The  main  feature  of  the  estimation  procedure  proposed  here  is  that  it 
is  recursive,  i.e.  as  the  radius  of  the  observation  disk  increases,  the  optimum 
estimation  filter  can  be  updated.  Earlier  methods  [9] ,  [10]  for  estimating 
isotropic  random  fields  were  nonrecursive.  An  advantage  of  this  procedure 
is  that  it  is  very  efficient:  it  requires  considerably  fewer  operations  (i.e. 
multiplications  and  additions)  than  methods  which  solve  directly  the  integral 
equation  associated  to  the  estimation  problem  over  a  finite  disk.  This 
property  is  not  surprising  if  we  note  that  the  vibrating  membrane  equations 
discussed  here  generalize  the  Levinson  recursions  [11],  [12]  and  the  vibrating 
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string  equations  of  [1]  -  [5]  which  are  known  to  provide  efficient  solutions 
of  the  cne-dimensional  estimation  problem  over  a  finite  interval. 

The  vibrating  membrane  equations  that  we  consider  are  identical  to 
equations  that  appear  in  the  Gelf and-Levitan  formulation  of  the  inverse 
scattering  problem  of  quantum  mechanics  [13]  ,  [14] .  By  transforming  the 
estimation  problem  over  a  finite  disk  into  an  equivalent  problem  in  the 
spectral  domain,  we  obtain  some  two-dimensional  Schrodinger  equations. 

These  equations  are  those  satisfied  by  partial  waves  of  nonzero  angular 
momentum  for  a  radially  symmetric  potential.  By  using  this  observation, 
we  show  that  only  one  vibrating  membrane  equation  needs  to  be  solved  to 
obtain  a  complete  solution  of  the  two-dimensional  estimation  problem. 

This  paper  is  organized  as  follows.  In  Section  II  we  expand  the 
observed  random  field  in  Fourier  series  and  introduce  some  one-dimensional 
filtering  problems  for  the  Fourier  coefficient  processes.  These  filtering 
problems  are  solved  by  exploiting  the  structure  of  the  covariance  kernels 
of  the  coefficient  processes,  and  in  Section  III  this  solution  is  used  to 
estimate  an  arbitrary  random  variable.  In  Section  III,  we  also  transform 
the  two  dimensional  estimation  problem  to  the  spectral  domain,  and  we  derive 
some  Schrodinger  equations.  In  Section  IV  a  set  of  special  transformations 
is  introduced  to  relate  the  solutions  of  these  Schrodinger  equations  and  to 
show  that  only  one  of  them  needs  to  be  computed.  The  construction  of  random 
field  estimates  is  considered  in  Section  V  and  the  special  case  when  we  want 
to  estimate  the  random  field  at  the  origin  of  the  observation  disk  is 
examined.  In  this  case,  the  optimum  filter  constitutes  an  approximation  of 
the  smoothing  filter  for  observations  over  the  whole  plane.  The  Section  VI 
presents  some  conclusions  and  discusses  some  possible  extensions  of  these 


results. 


1 


II .  Fourier  series  expansion  of  the  observed  field 

In  this  paper,  we  will  consider  the  estimation  problem  where  we  are 
given  some  observations 

dy(r,9)  =  z(r,6)  dA  +  dv(r.S)  (2.1) 

of  a  two-dimensional  isotropic  zero-mean  Gaussian  Random  field  z(*,*)  over 
a  disk  D  centered  at  the  origin  and  of  radius  R,  i.e.  for 

0  <  9  <  2m,  0  <  r  <  R  .  (2.2) 

Here  dA  =  rdrdQ  is  an  arbitrary  surface  element  and  v(*,*)  is  a  two-dimen¬ 
sional  Wiener  process  which  is  uncorrelated  with  z  ( • , • )  and  whose  incre¬ 
ments  are  such  that 

!dA  for  r=p  ,  9  =  6 

(2.3) 

0  otherwise 

Since  the  field  z(*,*)  is  isotropic,  its  covariance  is  given  by 

E[z(r,0)  z(p,d>)]  =  k(d)  (2.4) 

2  2 

where  d  =  (r  +  p  -  2rp  cos(9-<j>))  is  the  Euclidean  distance  between  the 
points  with  polar  coordinates  (r,6)  and  (p,4>).  If  z(*,*)  is  mean-square 
continuous,  or  equivalently  if  k(')  is  continuous,  k(*)  admits  a  spectral 
representation  of  the  form  (see  [15],  [16]) 

fee 

k(d)  =  i  J„ (Ad)  dM(X)  (2.5) 

Jo  G 

where  JQ(*)  is  the  zero-th  order  Bessel  function  and  where  M(*)  is  bounded 


and  nondecreasing. 
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The  Hilbert  spaces  spanned  by  the  observations  and  by  the  field  z(*,*) 

p 

will  be  denoted  respectively  as  y  '  =  K(dy(r,3) ;  o<r<R,  0<9<  2~)  and 

Z  =  H(z(r,S);  0<r<°°,  0<6<  2tt)  .  If  n  and  £  are  some  arbitrary  elements 

£ 

of  Y  and  Z,  they  can  be  represented  as 


dy(r,0) 


(2.6) 


£<r,9) 


z(r,8)  dA 


(2.7) 


where  r\(*,*)  and  £(•,*)  are  respectively  square-integrable  over  D  and  D 

R  00 

with  respect  to  the  measure  dA.  denotes  here  the  whole  plane,  or  equivalently 
a  disk  of  infinite  radius  centered  at  the  origin.  Note  that  while  the  observa¬ 
tions  are  given  over  D  ,  the  field  z(*,‘)  is  defined  over  D  . 

R  00 

The  estimation  problem  that  will  be  studied  here  consists  of  finding  the 
conditional  mean  of  a  random  variable  ^62  given  Y  .  A  special  case  that  will 
be  of  interest  is  when  £  =  z(Q)  is  the  field  at  the  origin.  The  optimum 
smoothing  filter  that  is  obtained  in  this  case  can  be  viewed  as  an  approximation 
for  the  optimum  smoothing  filter  based  on  observations  ever  the  whole  plane. 


A.  Fourier  coefficient  processes 

Since  the  field  z(r,9)  is  a  periodic  function  of  9 ,  it  can  be  expanded 
in  Fourier  series  as 

00 

z(r,9)  =  2  2  (r)  exp  jn9  (2.8) 

n=-co  n 


where  the  Fourier  coefficients 


z 

n 


(r) 


z (r,9) 


exp-jn0  d0 


(2.9) 


I 
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define  some  one-dimensional  stochastic  processes.  It  is  shown  in  [16]  that 
we  have 


E[z  (r)  z*(s)] 
n  m 


<*>  (r,s) 
n 


6 

nm 


wi  th 


5  (r,s)  =  |  J  (Ar)  J  (?\s)  d>t(A)  , 

n  l0  n  n 


(2.10) 


(2.11) 


where  J  (•/  is  the  n  '-order  Bessel  function  and  where  5  =  1  if  n  =  m 

n  nm 

and  6  =0  otherwise.  Thus,  the  crocesses  z  (•)  are  uncorrelated,  and  if 

nm  *  n 

we  denote  bv  2  =  H(z  (r) ;  o  <  r  <  °°)  the  Hilbert  space  scanned  bv  the  orocess 
'  n  n 

z  (•),  Z  can  be  decomposed  orthogonally  as 


00 

Z  =  Z  .  (2.12) 

n=-«>  n 


Since  the  covariance  of  the  processes  z  (•)  takes  the  form  (2.11)  it  is  clear 

that  z  (•)  is  not  stationary.  However,  the  covariance  d  has  the  following 
n  n 

property. 

Lemma  1.  Displacement  property  of  If  k(»)  is  twice  differentiable,  or 

equivalently  if  the  field  z(*,*)  is  mean-square  differentiable,  we  have 

-2  _  ~  2  .2  _  ^  2 

((-~  +  -  —  -  -  (— — r  +  -  —  -  —))  $(r,s)  =  0  (2.13) 

or  r  Sr  r  3s  s  9s  s 

with  the  boundary  conditions 


3r 


^0<r'  s) 

'  r=0 


0 


(2.14a) 


1 

*  * 
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and 


t>  (0,s)  =  0  for  n  ^  0  . 
'  n 


(2.14fc) 


Proof:  Since  k(*)  is  twice  differentiable,  6  (•,»)  is  also  twice  differentiable. 

-  n 

Then,  since  J  (Ar)  obeys  the  differential  eauation 
n 


••  1  *  9  n 

J  +-J  +  (X^  -  ~)  J  =0 
n  r  n  2  n 

r 


we  have 


32  .  1  3  n2 


__  + - _)  <*>  (r,s)  =  -  ; 

2  2  n  j 

r  r  3r  r  J  o 


r°°  „ 

!  Xz  J  (Ar)  J  (As)  dM ( A ) 
n  n 


~2  ,  *  2 
,3  13  n  ,  ,  ,  . 

=  ( — r  +  -  - - r)  (r,s) 

^  n 

9s  s  3s  s 

so  that  (2.13)  is  satisfied.  The  boundary  conditions  (2.14)  can  be  obtained 

bv  notina  that  J.(0)  =  0  and  that  J  (C)  =  0  for  n  /  0. 

'  u  n 


Thus,  although  z  (•)  is  not  stationary,  it  has  as  much  structure  as 
n 

a  stationary  process.  In  fact,  the  displacement  operator 

.2  ,  *  2 


6  =  {jL  +  I JL  _  si)  _  (11-  +  i  i_  _  Hi) 

n  3r2  r  3r  r2  3s2  s  3s  s2 


(2.15) 


associated  to  (2.13)  can  be  viewed  as  a  generalization  of  the  displacement 
operator 


2_  +  i_ 

3r  3s 


which  was  used  in  [11],  [12]  to  characterize  Toeplitz  kernels,  i.e.  covariances 
of  stationary  processes.  The  operator  6^  is  also  similar  to  the  one-dimensional 
wave  operator  introduced  in  [5] . 


The  noise  and  observation  processes  can  also  be  expanded  in  Fourier 
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so  that  the  processes  v^(*)  are  uncorrelated  Wiener  processes  cf  intensity 
r/2:r . 

A  consequence  of  this  property  is  that  the  estimation  problems  (2.18) 

£ 

can  be  solved  indeoendentlv  of  each  other.  Equivalently,  if  Y  =  H(dy  (r) ; 

'  n  n 

0  <  r  <  R)  is  the  Hilbert  space  spanned  by  the  n'h  observation  process,  by 

R  i  R 

combining  (2.10)  and  (2.19)  one  has  Y  _L  Y‘  for  n  ^  m.  To  show  that 

n  m 


yR  0  yR 

n=-oo  n 


(2.20) 


one  needs  to  prove  that  if 


n 


g ( r , ? )  dy ( r ,  3  ) 


(2.21) 


is  an  arbitary  element  of  Y  ,  then  n  can  be  expressed  as  the  sum  of  elements 
£ 

in  Y  for  n  =  0,  tl,  ±2  ...  But  p(r,8)  can  be  expanded  in  Fourier  series  as 
n 


oo 


n(r,0)  =  Z  n  (r)  exo  jn8 
n=-°°  n  * 

and  if  n(*,*)  is  real  g*(r)  »»  g  (r)  .  By  substituting  this  expression  inside 

(2.21)  and  by  usina  the  definition  (2.17)  of  v  (*)  one  finds  that 

n 


n  =  2f  Zm 

n=-“> 


0*(r)  dy  (r) 

h  n  n 


(2.22) 


so  that  g  is  the  sum  of  elements  of  Y  with  n  =  0,  ±1,  ±2...  This  proves  (2.20). 

n 


3.  Filtering  problem  for  the  Fourier  coefficients 


The  first  step  towards  the  solution  of  the  two-dimensional  estimation 
problem  described  above  is  to  find  the  filtered  estimate  of  the  field  z(R,9) 
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£ 

or.  the  edge  of  the  disk  D  given  the  observations  Y  .  Ey  using  the 

R 

R  R 

decomposition  (2.29)  of  Y‘  ,  the  conditional  mean  of  z(R,?)  given  Y  can  be 

expressed  as 


E[z(R,S) ! YR]  =  Z  E[z  (R) ,YR]  eXD  j 
'  n=-=°  n  n  - 


no 


(2.23) 


where 


f  P. 

E[z  (R) |y  ]  =  z  (r]r)  =  I  g  (R,s)  dy  (s) 
n  n  n  '  j  n  ■'n 

denotes  the  filtered  estimate  of  z  (R)  given  YR. 

n  n 

The  two-dimensional  filrering  problem  has  therefore  been  reduced  t 
set  of  one-dimensional  filtering  problems  for  the  Fourier  coef f icients . 
the  filtering  error  associated  to  (2.24)  is  denoted  as 


(2.24) 


o  a 


Tf 


Z  (RR)  =  Z  (R)  -  Z  (RiR)  , 
n  n  n 


(2.25) 


by  using  the  orthogonality  property  z^(r|r)  J-Yn  of  linear  least-squares 

estimates,  we  find  that  g  (*,*)  satisfies  the  integral  ecruatior. 

n 


fR  ! 

(R,r)  =  j  g  (R, s)  9  (s,r)  sds  +  —  g  (R,r) 
n  u  n  n  n 

'Q 


(2.26) 


for  0  <  r  <  R. 


Note  that  since  we  have  assumed  that  z(*,’)  is  mean-square  continuous. 


the  Fourier  coefficients  z  (*)  are  also  mean-square  continuous,  and  therefore 


$  (•,•)  is  continuous.  Thus,  the  operator 


fR 

<5  :  a(r)  -*■  b(r)  =  i  <i>  (r,s)  a(s)  sds 
-n  L  n 

•0 


is  defined  over  L  (rdr;  [0,R],  and  since  <J>  (*,*)  is  a  covariance  kernel,  the 


operator  is  self-adjoint  and  non-negative  definite,  so  that  +  I/2tv  is 
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invertible.  This  guarantees  the  existence  and  unicity  of  a  solution  in 

1*2  (rdr;  [0 , Ft ]  )  to  the  integral  equation  (2.26). 

To  compute  g  (*,*),  instead  of  solving  directly  the  integral  equation 

(2.26),  we  can  excloit  the  displacement  orooertv  of  C  .  This  oives  the 

■  n 

following  result. 


Theorem  1.  Vibrating  membrane  ecuation  for  a  .  If  the  kernel  k(*)  is  twice 
_ ~  n 

differentiable,  or  equivalently  if  z(*,*)  is  mean-square  differentiable,  gn ( * , * ) 
satisfies  the  partial  differential  equation 


+  -  —  -  \))  g  (R ,  r)  =  a  (R)  g  (R,r) 
^  2  .  2  -n  -n  'n 

or  r  or  r 


:?) 


with  the  boundary  conditions 


a  (R$ 
*n 


“  dR 


(Ra  (R,  R)  ) 
n 


(2.28) 


and 


g  (R,r)|  =0,  g  (R,0)  =  9  for  n  ^  0  .  (2.29) 

or  o  _  n 

r=0 

The  proof  of  Theorem  1  is  given  in  Appendix  A.  The  equation  (2.27)  can 
be  viewed  as  being  satisifed  by  a  vibrating  membrane  submitted  to  a  circularly 
symmetric  excitation.  For  such  a  membrane  R  and  r  denote  respectively  the 
radial  position  and  time.  Since  (2.27)  takes  nearly  the  same  form  as  the 
displacement  relation  satisfied  by  ®  ,  it  can  be  viewed  as  a  generalisation 
of  the  Levinson  recursions  [11] ,  [12]  and  of  the  vibrating  string  equations 

introduced  in  [5]  for  the  one-dimensional  estimation  problem.  Note  indeed 
that  these  recursions  were  obtained  by  exploiting  the  displacement  properties 
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of  Toeplitz,  and  of  Toeplitz  plus  Hankel  kernels. 

By  performing  the  transformation 
1/2 

l  (R,r)  =  (Rr)  7  g  (P.,r)  (2.30) 

n  n 

the  recursions  (2.27)  can  be  expressed  as 
2  2 

( — r  -  -4  (n2  -  f))  -  fi—r  -  4r  (n2  -  \))l  (R,r>  =  o  (R)  Z  (R,r)  (2.31) 
^  2  2  4  »  2  2  4n  *n  n 

or  R  3r  r 


with  the  boundary  condition 


c  (R) 

■■•ri 


o  A. 

“  dR 


i  (R,R)  . 

n 


This  equation  appears  in  the  study  of  inverse  scattering  problems  of  quantum 

mechanics  [13] ,  [14]  and  the  relation  between  these  problems  and  the  linear 

estimation  results  discussed  here  will  be  explored  in  [17].  The  filter  Z  (*,*) 

n 

can  be  interpreted  by  considering  the  normalized  estimation  problem 


dv  (r)  =  z  (r)  dr  +  dv  (r)  ,  0  <  r  <  R  (2.32) 

n  n  n  - 

-1/2  -  -1/2  -  1/2 
where  dy  (r)  =  r  dv  (r)  ,  dv  (r)  =  r  /  dv  (r)  and  z  (r)  =  r  '  z  (r)  . 
n  *  n  n  n  n  n 

In  this  case,  the  noise  process  v  (•)  has  stationary  increments,  and  ^(R,*) 

~R  — 

is  the  optimal  filter  for  estimating  z^CR)  given  Yn  =  H(dyn(r) ,  0  <  r  <  R) . 

The  mean-square  error  associated  to  the  filtering  problem  (2.18)  can  be 
denoted  as 

e2(R)  =  E[z2(R|R) ] 
n  n  1 

and  by  using  the  orthogonality  of  linear  least-squares  estiamtes  we  find  that 

e2  (R)  =  g  (F.,R)  . 
n  n 


(2.33) 
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This  relation  provides  an  interpretation  of  the  potential  q^(R)  given  by 

(2. 28)  . 

C .  Innovations  processes  for  the  Fourier  coefficients 


Given  the  observations  (2.18),  we  have  shown  how  to  compute  the  filtered 

estimates  z  (r,r).  These  estimates  can  be  used  to  aenerate  an  innovations 
n 

process 

dv  (r)  =  dy  (r)  -  z  <r|r)  dr  =  z  (r|r)  dr  +  dv  (r)  (2.34) 

n  -'n  n  n  n 


which  is  a  Wiener  process  of  intensity  r/2i ,  i.e. 


E  [dv  (r)  dv*  (s)  ] 
n  n 


~  dr  for  r=s 

0  otherwise 


(2.35) 


Furthermore,  the  Hilbert  space  v  =  H(dv  (r)  ,  0  <  r  <  R)  generated  by  this 

n  n  - 

R  til 

process  is  identical  to  the  Hilbert  space  Y^  associated  to  the  n  Fourier 
coefficient  of  the  observations,  so  that  the  Hilbert  space  decomposition 
(2.20)  reduces  to 


Y 


R 


00 


© 


(2.36) 


III .  The  general  estimation  problem 

We  have  shown  above  how  to  solve  the  filtering  problem  associated  to 
the  Fourier  coefficients  of  the  observed  field  z(*,*).  This  result  will  now 
be  used  to  estimate  an  arbitrary  random  variable. 


Theorem  2.  Let  £  be  a  zero-mean  random  variable  whose  joint  statistics  with 
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variables  in  Y  are  known  and  are  Gaussian,  A  special  case  is  cf  course 

p 


when  •  S  I.  Then,  the  conditional  mean  and  error  variance  of  %  given  V" 
are  civen  bv 


.R, 


30  ,'R 


•  19  !  l 


n=-°%  n 


(3.1) 


and 


=  E[(E  -  E  [?|yR])2]  =  E[^]  -  2tt  J  m  I  |E  (r)|*  rdr 


fR,_  .  .,2 


n=-“  J0  1 


(3.2) 


where  £,  (•)  is  such  that 


E[£  d'J*  (r)  ]  =  C  (r)  rdr  . 
n  n 


(3.3) 


Proof:  Use  the  decomposition  (2.36)  and  note  that  v  (•)  is  a  Wiener  orocess  of 
-  n 


intensity  r/2rr. 


The  previous  results  can  be  expressed  entirely  in  the  spectral  domain. 

To  do  so,  we  denote  by  Y  =  H(dy(r,9);  0  <  r  <  °°,  0  <  6  <  2tt)  the  Hilbert  space 
spanned  by  the  observations  over  the  whole  plane,  and  we  use  the  Kolmogorov 
isometrv 


f 


_  f 


a  =  |  a(r,9)  dy(r,9)  a (X , d) )  =  !  a(r,0)  exp(jXr  cos(0-<M)  rdrd9  (3.4) 


Deo 


between  Y  and  L„(dF)  where  F(X,<J>)  is  the  spectral  measure  associated  to  the 
observations  y(*,*).  This  measure  is  given  by 


k  (r)  + 


6  (r) 


irr 


exptjXr  cos(9-4>))  dF(X,>$) 


(3.5) 


where  <5  ( r)/rr  is  a  two-dimensional  impulse  function  written  in  polar  coordinates. 
Note  that  in  (3.4)  and  (3.5)  the  term  Xr  cos(9-d>)  corresponds  to  the  inner  product 


i 
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betveen  the  vectors  x  =  (r  cos  6,  r  sin  9)  and  u  =  (X  cos  <J> ,  X  sin  >■)  . 

Since  the  left  hand  side  of  (3.5)  depends  on  r  only,  dF ( • , • )  is  a 
function  of  X  only,  i.e. 


dF  (X  ,  £ )  =  dF  (X) 


d<f> 


2m 


so  that  (3.5)  reduces  to 


(3.6) 


k(r> 


6(r) 

mr 


Jn(Xr)  dF  (X) 
Jo  0  1 


(3.7) 


where  we  have  used  the  fact  that 


J  (Xr) 


2TT 


J-2TT 

!  exp(jXr  cos  (8—i>) )  d$ 
JQ 


for  6  arbitrary.  Furthermore,  F^(')  can  be  expressed  in  function  of  the 
measure  M (• )  introduced  in  (2.5)  as 


dF1(X)  =  dM(X) 


XdX 

2m 


(3.8) 


The  isometry  (3.4)  is  such  that  if  a  and  b  are  two  random  variables 
of  Y  corresponding  to  some  functions  a(X,<J)  and  6(X,i)  of  L  (dF)  ,  one  has 


Stab]  =<a,S>p=  a  (X  ,<i> )  £*(X,$)  dF(Xr<p ) 


(3.9) 


In  this  framework,  the  subspace  Y  of  Y  generated  by  the  observations  over  D 

P 

£ 

is  isomorphic  to  the  subspace  S  of  L2(dF)  generated  by  the  functions 


a(X,<>)  =  )  a(r,9)  exp(jXr  cos  (8-<j>))  rdrdQ 

D_ 


(3.10) 


where  a<  *,  *)  is  square-integrable  over  D  with  respect  to  rdrd6. 

K 


i 

J 


A.  Decomposition  of  S 


To  obtain  a  decomposition  of  S  equivalent  to  the  decomposition  (2.20) 
for  Y  ,  we  note  chat  if 


n  =  —  j  n(r)  exo-jn0  dy(r,9) 
“  -dr 


(3.11) 


Ls  an  arbitrary  random  variable  of  Y‘  ,  h  is  mapped  into  the  function 

n 


ft(A,£)  =  ~  [  h(r)  exp  j  (Xr  cos  ( 9— )  -  n9)  rdrdQ . 

J°R 


(3.12) 


By  taking  into  account  the  Fourier  series  expansion 


exp  (jx  cos  9)  *  L  j*  J  (x)  exo  jn9 
n=-°°  n  * 


^(•,*)  can  be  expressed 


ft  (A  ,$)  =  ft1(X)  exp  -  j ng> 


(3.13a) 


where 


ft,  (X)  =  j  h(r)  J  (Xr)  rdr 

1  Jo 


(3.13b) 


R  R 

so  that  is  isomorphic  to  the  subspace  Sn  of  (3F)  generated  by  functions 

R  i  R 

of  the  type  (3.13).  From  (3.13a)  it  is  clear  that  S  JL  S  for  n  ^  m.  Moreover, 

n  m 

/V  £ 

if  ft(*,*)  and  £(*,*)  are  two  functions  of  S  ,  we  have 

n 


<  ft.  t  >„  =  <  n.,  5.  >  =  n,  (X)  K*  (X)  dp.  (X)  , 

1  1  Fi  J0  1  1  1 


(3.14) 
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and  the  mapping 


M\,})  =  t  (\)  exp  -  jn6  — ►  h,  (X )  (3.15) 

R  R 

is  an  isometrv  between  S  and  the  subsoace  L'  of  L_ (d?  )  generated  by  functions 

n  n  2  x 

of  the  form  (3.13b). 

3y  transforming  the  decomposition  (2.20)  under  the  isometry  (3.4),  we 
find  that 

00 

SR  =  ©  SR  .  (3.16) 

n 

n — oo 


In  the  special  case  when  R  =  °°,  the  decomposition  (3.16)  becomes 


oo 


(3.17) 


where  S  is  the  subsoace  of  L„(dF)  sDanned  bv  the  functions  h,  (A)  exp  -  ini', 
n  2  1 

/v 

where  h^(’)  belongs  to  L2(dF^).  Note  that  if  h^(*)  is  arbitrary  function  of 
L2^<iRl''  ^  Can  exPresse<*  as 


r, 

1 


(A) 


h  (r) 


J  (Ar) 


rdr 


where  h ( • )  is  the  nth  -order  Hankel  transform  [18] 


of  n1  ( • )  . 


B.  Orthogonalization  of  the  functions  (Ar) 


A.  p 

The  relation  (3.13b)  shows  that  if  n, (•)  belongs  to  L  ,  it  can  be 

1  n 

expressed  as  linear  combination  of  the  functions  (j  (Xr)  ,  0  <  r  <  r)  . 


* 
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However,  since  we  have 


<  J  (Ar),  J  (As)  >_  =  $  (r-s)  +  9  (r,s) 

n  n  F  2"rr  n 


(3.13) 


these  functions  are  not  orthogonal  in  L^tdF^).  To  orthogcnalize  then,  we  will 
use  the  correspondence 


rR 

Jo 


a  (r)  dv  (r) 
n 


.n 


a  (r)  Y  (r,A)  rdr 
n 


(3.19) 


R  R 

between  Y  and  L  ,  where 
n  n 


Y  (r , A)  =  J  (Ar)  -  f 

n  n  Jo 


g  (r,s)  J  (As)  sds 
n  n 


(3.20) 


and  where  a(*)  is  an  arbitrary  function  of  (rdr;  [ 0 , R] ) .  Then,  since  the 

« 

innovations  process  V  (•)  is  a  Wiener  process  of  intensitv  r/2T,  the  functions 

n 

Y  (• ,A)  are  orthogonal  and 


c  yn(r,A),  yn(s,A)  ^  5  (r-s)  . 


(3.21) 


To  characterize  the  functions  y  (*,*),  instead  of  using  the  functions  g  (*,•), 

n  '  n 

we  can  use  the  following  result. 


Theorem  3.  For  n  >  0,  y  (*,A)  satisfies  the  differential  equation 
-  n 


2  ,  2 

(_1_  +  -  —  +  (A2  -  q  (r)  -  ^r-))  y  (r,A)  =  0 
«  n  £  n 

dr  r  dr  r 


(3.22) 


with  the  boundary  condition 


lim  2n  n!  (Ar)  n  y  (r,A)  =  1; 

_  n 

r+0 


(3.23) 


-20- 


ar.d  for  n  <  0,  we  have 


Y  (r,\) 
n 


(-1)  Y_n<r,\) 


(3.24) 


Proof :  By  operating  with 


...  d  la  n 

A  (r)  =  — r  + - 7 

-  n  .2  ,  2 

dr  r  ar  r 


(3.25) 


on  (3.20),  one  gets 


A  (r)  y  (r,A)  =  -  (A ^  -  a  (r) )  J  (Ar)  +  r  ■£—  a  (r,s)  i  J  (Ar) 
-n  n  n  n  os  'n  n 


dJ  [r 

-  Xrg  (r,r)  — (Ar)  -  j  (A  (r)  g  (r,s))  J  (As)  sds 
n  dr  J  -n  Tn  n 


Then,  since  the  displacement  operator  can  be  exoressed  as  £  =  A  (r)  -  A  (s)  , 

-n  *  -n  -n  -n 

by  using  the  aquation  (2.27)  for  gn(r,s)  and  integrating  by  parts,  we  obtain 

(3.22).  The  boundary  conditions  (3.23)  and  (3.24)  are  identical  to  those 

satisfied  by  J  (Ar)  and  can  be  derived  from  (3.20).  ■ 

n 

The  differential  equation  (3.22)  takes  a  familiar  form  if  we  perform 
the  transformation 


7  (r , A)  =  (rA)1/2  y  (r,A) 
n  n 


(3.26) 


In  this  case  (3.22)  becomes 


7  (r, A)  +  (A2  -  q  (r)  -  ~  (n2  -  -7) )  7  (r,A)  =  0 
n  n  4  n 

r 


(3.27) 


with  the  boundary  condition 


lim  2  n!  (Ar) 
r-K) 


7n(r,A)  =  1 


- (n+1/2) 
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for  n  >  0,  and  (r,.\)  =  (-l)n  (r,A)  for  n  <  0.  For  a  space  of  dimension  two, 

n  -n  ~ 

this  is  the  radial  Schrodinger  equation  satisfied  by  the  partial  wave  of 

angular  momentum  r.  associated  to  a  particle  moving  in  the  circularly  symmetric 

ootential  a  (•)  (see  [17]  for  more  details). 

'  *n 

C .  Projection  on  S 

Under  the  isomorphism  (3.4),  if  5  is  an  arbitrary  random  variable  of  Y 
corresponding  to  the  function  £(A,$)  of  L  (dF)  ,  the  problem  of  finding  the 
conditional  mean  of  %  given  Y*  is  equivalent  to  the  one  of  projecting  ?(A,i>) 
on  3  .  Thus,  we  have  the  correspondence 


E(g|YR]  P  5U,$) 

S 

where  P  denotes  the  oroiection  operator  on  S  .  Then,  if 
SR 

CO 

£(A,$)  =  Z  c  (A)  exp  jn<{> 

1-1=-°°  n 


(3.28) 


(3.29) 


denotes  the  Fourier  series  expansion  of  £(•,•),  Theorem  2  can  be  reformulated 
as  follows. 

^  ^  R 

Theorem  2  ;  The  projection  of  £(A,4>)  on  S'  is  given  by 


P  8  €  (A,4» 


5  (r) 

n 


Y  (r,A)rdr 
n 


’) 


exp 


jnd 


(3.30) 


where 


^n(r) 


c  ~ 

5  (X)  Y*  (r ,A)  dF. (A) 

-n  n  1 

0 


(3.31) 


satisfies  (3.3). 


Proof:  By  replacing  n  by  -n  in  (3.29),  the  expansion  (3.29)  can  be  viewed 

A 

as  decomoosition  of  ?(•/)  in  terms  of  elements  of  S  with  n  =  0,  ±1,  i2... 

n 

R 

The  decomposition  (3.16)  for  S  can  then  be  used  to  show  that 

00  ^ 

P  _  C(A,J)  =  -  P  „  (  CD  sxd  -  jnti 

R  n=-°°  R  -n 

S  Li 

n 

where  P^p  denotes  the  projection  operator  on  the  subspace  L?-  of  L^tdF^). 

n  11  R 

Furthermore,  by  using  Theorem  2  and  the  isometry  (3.19)  between  Yp-  and  L^', 

we  find  that 


/\ 


R 


•R 

2tj"  |  ?  (r)  y  (r,,\)  rar 

J0  n  n 


where  ?n D )  is  given  by  (3.3).  The  relation  (3.31)  is  equivalent  to  (3.3) 
under  the  isometry  (3.19).  B 


Note  that  the  solution  of  the  general  estimation  problem  that  we  have 

obtained  in  Theorems  2  and  2'  deoends  on  the  functions  (g  ,  n  6  J  !  or  {'(  , n  C  v)  . 

n  n 

It  will  now  be  shown  that  it  is  not  necessary  to  compute  all  these  functions,  and  that 

for  |n|  >  1,  g  and  y  can  be  expressed  in  function  of  gn  and  yn. 

“  n  n  u  u 

IV.  Special  transformations 


The  first  step  in  order  to  relate  g  and  a  (respectively  y  and  y  )  is 

n+i  "n  n+J.  n 

to  note  that  the  covariance  kernels  d>  (*,*)  satisfy  the  followina  property. 

n 


Lemma  2.  If  k(«)  is  differentiable,  for  all  n  we  have 


(£  -  b  K+i{x’a)  5  ° 


\  /  Y  _ 

3r  r  n 


(—--”)<(»  (r , s)  +  (vr  +  --n+-1-)  4>_i  (r,s>  =  o 
3s  s  n  or 


n+1 


(4.1a) 


(4.1b) 


r 
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Proof  :  The  Bessel  functions  iJ  (r) ,  n  6  Z!  satisfy  the  recursions 


J  (r)  +  -  J  (r)  =  J  .  (r) 
n  r  n  n-1 


(r)  -  -  J  (r)  =  J  (r) 
i  r  n  n+1 


and  A  (•,*)  is  given  bv  (2.11).  This  imolies  that 
n 


(4-  -  -)  *  (r ,  s)  = 

a  r  r  n 


Xj  , ,  (Xr)  J  (Xs)  dM  (X ) 


Jo  n+1 


and 


(^-  -  (~— )  <P  .  (r, s)  =  i  X J  xl(Xr)  J  (Xs)  dM(X) 
:s  s  n+1  jQ  n+1  n 


so  that  (4.1a)  is  satisfied.  The  relation  (4.16)  can  be  obtained  by  symmetry. 


Note  that  if  A  (r)  is  civen  bv  (3.25),  we  have 

-n 


V  ■  ^+i2r1'  &-§> 


(4.2a) 


,  ,  \  ,9  n,  ,  3  (n+1) , 

A  . .  (r)  =  (- - )  (r—  +  - ) 

-n+1  or  r  dr  r 


(4.2b) 


and  since  6  =  A  (r)  -  A  (s) ,  the  displacement  propertv  (2.13)  for  i  can  be 

-n  -n  -n  *  n 

obtained  by  combining  (4.1a)  and  (4.1b).  Then,  by  using  the  integral  equation 
(2.26)  for  gn»  we  obtain  the  following  result. 


Theorem  4.  If  k(’)  is  differentiable,  the  functions  g  (• 
-  n 


r  n 


3  (n+1), 

'3s 


-  3  g^(r,s)  +  (fz  +  g_ , ,  (r,s)  = 


'n+1 


-  a  (r) 
n 


,  • )  are  such  that 

g  (r,s)  (4.3a) 

n 
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— )  g  (r,s) 
s  n 


+ 


(n+l)} 

r 


a  , ,  (r ,  s)  =  a  (r)  a  ,  (r , s) 
•n+1  n  'n+1 


(4.3b) 


wr.ere 


a  (r)  =  r  (g  (r,r)  -  g  (r,r)).  (4.4) 

n  n  n+1 

The  proof  of  Theorem  4  is  giver,  in  Appendix  A.  This  result  can  be  used 
to  derive  the  vibrating  membrane  equation  (2.27),  since  by  combining  (4.3a)  and 
(4.3b)  and  using  the  identities  (4.2),  we  obtain  (2.27)  with 


,  .  2 ,  ,  .  ,  .  (2n+l)  ,  , 

q  (r)  =  a  (r)  -  a  (r)  -  -  a  (r) 

n  n  n  r  n 

7  (2n-*-l) 

a  (r)  =  a  (r)  +  a  (r)  -  -  ■  a  (r)  . 

n+1  n  n  r  n 


(4.5a) 

(4.5b) 


The  relations  (4.5)  show  that  the  potentials  a  (•)  and  a  , (•)  are  related, 

-n  -r.+l 

since  they  can  be  expressed  in  terms  of  the  same  function  a  (*).  Converselv, 

n 

civen  a  (•)  or  a  , (•)  we  can  obtain  a  (•)  bv  solving  the  Riccati  eouation  (4.5a) 
*n  'n+1  n 

or  (4.5b)  with  the  boundary  condition 

a  (0)  =  0  for  all  n. 
n 


The  functions  a  (•)  play  here  a  role  similar  to  the  reflection  coefficient  function 
n 

of  [5]  . 

An  alternative  method  of  computing  a  (•)  aiven  either  a  (•)  or  c  , (•)  is  to 

n  -n  n+1 

substitute 


a 

n 


(r) 


n 

r 


w  (r) 
n 

w  (r) 
n 


or 


a  (r) 
n 


(n+1) 

r 


u  (r) 

+  -a - 

u  (r) 
n 


(4.6a) 


(4.6b) 
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wit  h  the  boundarv  conditions 


lim  2n  n!  r“n  w  (r)  =  1 

n 

r-K) 


(4.7a) 


yj-l 

lim  r  u  (r)  =  1 

»— *O0  ^ 


(4.7b) 


inside  (4.5a)  or  (4.5b).  Then,  the  equations  (4.5a)  and  (4.5b)  are  linearized 
and  take  the  form 


2  2 
,  d  Id  ,  ,  .  n  .  .  .  .  . 

( — q  + - (c  (r)  +  — )  )  w  (r)  =  0 

.  2  n  2  n 

dr  r  dr  r 


(4.8a) 


,d  .Id  ,  ,  <  (n+1)  ,  , 

( — 0  - (q  ,  i  (r)  +  - g—  ) )  u  (r)  =  0 

.  -  .  n+1  l  n 

or  r  dr  r 


(4.8b) 


This  shows  that  w  (r)  and  u  , (r)  are  two  solutions  of  the  differential  ecuation 
n  n-1 

(4.8a).  However  w^(r)  is  regular  at  the  origin,  whereas  by  subtracting  (4.6a) 
from  (4.5b)  and  integrating,  we  find  that 


u  (r)  =  - f-r- 

n  rw  (r) 


(4.9) 


where  is  a  nonzero  constant.  By  taking  into  account  (4.7a)  we  see  therefore 
that  u^(r)  has  a  singularity  of  order  n+1  at  the  origin.  The  identity  (4.9) 
also  shows  that  the  singular  solution  associated  to  the  potential  q  ^(*)  can  be 
obtained  from  the  regular  solution  associated  to  q^(*).  We  can  identify 


w  (r)  =  9  (r, 0) 
n  n 


(4.10) 


where  9n(r,.\)  =  Yn(r,.\)/X  ,  and  by  integrating  (4.6a)  we  find  that 


w  (r)  =  - 

n  2°  n! 


fr 

exp  -  |  a  (u)  du 


Jo1" 


(4.11) 


6- 


The  identities  (4.3)  of  Theorem  4  can  be  rewritten  as 


r“  ~  (r_r‘  c  (r,s)  ) 
r  'n 


■  n4-!)  4  ,  n-rl  .  .  ,  .  ,  , 

—  (s  c  (r ,  s)  )  =  -  a  (r  c  (r,s) 
?s  'n+1  r.  n 


r  j  —  n 

s"  - —  (s  g  (r, s) )  + 
os  n 


-(n+x)  .•  ,  n-1  ,  .  ,  ,  , 

r  —  (r  o  (r,s) )  =  a  (r)  g  (r,s) 

or  "n+1  n  n+1 


and  therefore,  by  integration,  gn+-|  can  be  expressed  in  function  of  g  ana 

vice-versa.  However,  a  simpler  set  of  identities  can  be  obtained  by  considering 

the  functions  '(  (*,*). 

n 


rheorem  4  .  For  all  n,  we  have 


AY  ,(r,x>  =  (^  -  7)  Y  (r,\)  -  a  (r)  y  (r,X) 
n-ri  ar  r  n  n  n 


(4.12a) 


(n-4-!) 


■'■Vn(r,X)  =  (—+—7—)  Yr^(r,\)  -  an(r)  7n+J  (r,\) 


'  n-1 
d  n 


n  ’n+1 


Proof :  By  operating  with  ^  ^  on  (3.20),  one  gets 


(4.12b) 


(-7-  -  — )  Y  (r,X)  =  XJ  .  (Xr)  -  rg  (r,r)  J  (Xr) 
dr  r  n  n+1  n  n 


(4.13) 


(■  r 

(m—  -  — )  a  (r,s)  J  (Xs)  sds 
J0  5r  r  -n  n 


Then,  if  we  use  (4.3a)  and  integrate  by  parts,  we  find  that 


(7 - — )  g  (r , s)  J  (as)  sds  =  -  g  (r,s)  J  (Xs)s 

3r  r  n  n  n+1  n 


I  s=r 


t  s=0 


f3 


f3 

(r)  ! 


X  i  g  (r,s)  J  ,  (Xs)  sds  -  a  (r)  g  (r,s)  J  (Xs)  sds 
J  o  n+1  n+1  n  J  Q  n  n 


By  substituting  this  expression  inside  (4.13),  and  taking  into  account  (4.4) 
and  (3.20),  we  obtain  (4.12a).  The  identity  (4.12b)  can  be  proved  similarly. 
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When  tha  relations  (4.6)  for  w  (•)  and  u  (•)  are  taken  into  account, 

n  n 

the  identities  (4.12)  can  be  expressed  as 


Ay  (r,A)  =  WE’,'  ( r , .'v )  ,  w  (r)]/w  (r) 
n+l  n  n  n 


(4.14a) 


Ay  ( r , A )  =  W[y  (r,A) ,  u  (r)]/u  (r) 
n  n+l  n  n 


(4.14b) 


where 


W[f(r),  g(r)]  A  f (r)  g(r)  -  f (r)  or (r) 


denotes  the  Wronskian  of  f(*)  and  g(*).  This  class  of  transformations  was 
introduced  by  Crum  [19]  and  was  latar  used  by  Krein  [20]  (see  also  [13],  [14]) 
to  study  the  inverse  scattering  problem  of  quantum  mechanics  for  waves  of 
nonzero  angular  momentum.  Under  such  transformations,  the  solutions  of  the 
Schrodinger  equation  for  a  potential  q(»)  and  angular  momentum  n  can  be  trans¬ 
formed  into  solutions  for  a  potential  q'(*)  and  angular  momentum  n' ,  where  q'(*) 
behaves  like  q(*)  near  r  =  0  and  as  r  -*■  °°.  For  the  problem  consider  here,  we 
can  use  these  transformations  to  express  Y  (*,*)  in  function  of  Yg(*,*)  for  all  n 

However,  the  relation  (4.14a)  alone  is  not  sufficient  to  show  that  Y  (*,*) 

n 

can  be  expressed  in  function  of  Yq(»,*)  for  all  n.  Since  (4.14a)  depends  also 

on  w  ( • ) ,  we  need  to  prove  that  w  , , ( • )  can  be  computed  in  function  of  w  ( • ) . 
n  n+l 

To  do  so,  note  that  the  Wronskian 

W ( r , A )  =  W[y  <r,A)  ,  w  (r)  ] 
n  n 

satisfies  the  differential  equation 

4-  „  +  ~  W  =  -A2  y  (r, A )  w  (r)  (4.15) 

dr  r  n  n 
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with  W(0,A)  =  0  for  all  n.  3v  integration,  one  obtains 


( r , A )  =  -  ~ —  v  (s,.\)  v  (s)  sds 

r  ;0  n 

and  by  using  (4.14a),  we  find  that 


wn+l(r)  =  lin|  Vl(rA)/A 

A-^-0 


n-rl 


rw  (r) 
n 


w  (s)  sds 


(4.16) 


This  shows  that  y  (•,*)  is  a  function  of  y.i*,*),  so  that  we  need  to  solve 
n  0 

only  one  of  the  Fourier  coefficient  filtering  problems  associated  to  our 
original  two-dimensional  estimation  nroblem. 


V.  Random  field  estimates 


In  Section  III,  it  was  shown  how  to  estimate  an  arbitrary  random  variable 

i  whose  joint  statistics  with  y(*,')  are  known.  Vie  will  now  consider  the  case 

when  g  =  z(r,6).  In  this  case,  instead  of  using  the  estimation  method  described 

R 

in  Section  III,  we  can  decompose  the  conditional  mean  of  z(r,9)  given  Y  as 


E [z  (r ,6 )  ; YR]  =  -  E[z  (r) | YR]  exp  jn0 

n=-co  n  1  n 


(5.1) 


and  then  denote 


z  (r |R)  =  E [z  (r) |YR] 
n  n  n 


H  (r,s;  R)  dy  (s) 
n  n 


(5.2) 


The  problem  of  estimating  z(r,8)  given  Y  is  now  reduced  to  the  one  of 

computing  H  (•,  • ;  R)  for  all  n. 
n 

By  using  the  orthogonality  property  of  linear  least-squares  estimates, 
we  find  that  •;  R)  satisfies  the  integral  equation 


$  (r, s) 
n 


H  (r,u;  R)  $  (u,s)  udu  +  —  H  (r,s;  R) 
n  n  2ir  n 


(5.3) 


A 


A 


_ ( 


-29 


where  0  <  s  <  R  and  0  <  r  <  °°.  In  the  special  case  when  0  <  r  <  R,  i.e.  when 

the  field  :(•,•)  is  estimated  at  a  point  located  inside  the  disk  D  .  if  H  is 

R  -n 

the  operator  cn  L_,  (rdr;  [0,R])  defined  by 

•  p, 

H  :  a(r)  -  b(r)  =  '  H  (r,s;  R)  a(s)  sds 
J  0  n 

the  equation  (5.3)  car.  be  rewritten  in  operator  notation  as 

(1/2"  +  5  )  (I  -  H  )  =  1/2"  .  (5.4) 

-n  -  -n 

This  shows  that  H  is  the  resolvent  operator  associated  to  $  .  Furthermore, 

-n  1  -n 

by  setting  r  =  R  in  (5.3),  we  can  identify 


H  ( ?. ,  s ;  R)  =  g  (R,s)  .  (5.5) 

n  '  n 

To  compute  the  kernel  H  (*,  •;  R) ,  instead  of  solvincr  directlv  the  integral 

n 

equation  (5.3),  we  can  use  the  3ellman-Krein  identity  [21] 


1  1_ 

R  jR 


H  ( r ,  s  ; 
n 


R)  =  -  H  ( 
n 


,  R;  R)  a  (R,s) 
n 


(5.6) 


which  is  obtained  by  taking  the  partial  derivative  of  (5.3)  with  respect  to  R 

and  by  comparing  the  resulting  equation  to  (2.26).  By  solving  first  the  vibrating 

membrane  equation  (2.27)  for  g  (*,•)  we  can  use  this  identity  to  comcute  H  (•,  •;  R) 

n  n 

for  increasing  values  of  R.  In  the  special  case  when  0  <  r  <  R,  the  identity 
(5.6)  becomes 


i^Hn(r,s;  R)  -  -  gn<R,r)  gn(R,s>  (5.7) 

and  by  integration,  it  can  be  expressed  in  operator  form  as 

(I  -  H  )  =  (I  -  o  *)  (I  -  g  )  (5.8) 

-n  -n  -  ‘n 

where  g  is  the  lower  trianaular  Volterra  operator  on  L  (rdr;  [0,R] )  given  by 

'  *■ 

fr 

g  ••  a(r)  ■+■  i  g  (r,s)  a(s)  sds 
*n  n 

J0 


*  v>-  xa  -  w— •>  mr-e *■ *  t 
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and  where  a*  is  the  dual  operator  of  g  ,  i.e.  its  kernel  is  given  by 
-n  £n 

c*  (r , s)  =  g  (s,r)  (5.9) 

n  r. 

for  all  C  <  r,  s  <  R.  The  operator  factorization  (5.3)  in  terms  of  upper  times 
lower  triangular  operators  is  of  the  type  considered  by  Gohberg  and  Krein  [22] . 


Comments 


1)  To  implement  numerically  the  previous  method  for  computing  H  (*,  *>  P) , 

the  interval  [0,R]  can  be  divided  into  N  subintervals  of  length  A  =  r/?j  and  the 

Bellman-Krein  identity  (5.7)  and  the  vibrating  membrane  equation  (2.27)  can  be 

discretized.  The  resulting  recursions  enable  us  to  ccmoute  H  (i A,  j A;  R)  for 

n 

all  0  <  i,  j  <  N  with  only  0(M3)  operations,  whereas  by  discretizing  the  integral 

equation  (5.3)  and  by  solving  the  corresponding  system  of  linear  ecuations,  we 
4 

would  need  0(N  )  operations.  This  shows  that  the  vibrating  membrane  ecuations 
of  Section  TI  provide  an  efficient  solution  of  the  estimation  problem  over  a 
finite  disk. 

R 

2)  The  expression  (5.1)  -  (5.2)  for  the  estimate  of  z(r,9)  given  V  is 
different  from  the  one  that  was  obtained  in  Theorem  2.  To  relate  these  two 


[  ! 


expressions,  note  that  the  integrated  form  of  the  Bellman-Krein  identity  (5.6) 


H  (r,s;  R)  =  H  (r,  s;  s)  -  \  H  (r,  u?  u)  g  (u,s)  udu 
n  n  in  n 

's 


(5.10) 


By  substituting  this  relation  inside  (5.2) ,  we  find  that 


z  (riR)  =  i  H  (r,  s;  s)  dv  (s) 

n  J0  n  n 


(5.11) 


so  that 


E[z(r,9)  dv*(s)J  =  ~  H  (r,  s*  s)  exp  jn0  sds  , 
n  n 


(5.12) 
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and  by  taking  (5.12)  into  account  inside  (4.1)^  and  (5.11),  we  obtain  the 
expression  (3.1)  for  the  estimate  of  5  =  z(r,®). 

The  symmetric  smoothing  problem 

A  special  case  of  interest  is  when  we  want  to  estimate  the  random  field 

z(0)  at  the  orioin  of  the  disk  D  .  In  this  case,  the  aecmetrv  of  estimation 
-  R  ' 

is  circularly  symmetric,  and  it  will  be  shown  below  that  the  associated  smoothing 

filter  can  be  used  to  approximate  the  smoothing  filter  based  on  observations 

over  the  whole  olane.  For  this  problem,  we  have  z(0)  =  zo(0)  and  z  (0)  =0 

for  n  ^  0,  so  that  the  estimate  of  z(0)  given  V  needs  only  to  be  based  on 
R  . 

Yq,  i.e. 

E  [z  (0)  j  YR]  =  zQ  (0  j  R)  .  (5.13) 

3v  denoting 

f  R 

z0 (0 i R)  =  j  f(R,r)  dy  (r)  (5.14) 

J  Jo  L 

=  f  (R,r)  dy  (r,9) 

dr 

we  find  that 


f (R,r)  =  Hq(0,  r;  R)  (5.15) 

and  by  using  (5.6)  one  gets 

|  f (R,r)  =  -  f (R,R)  g0(R,r)  (5.16) 

with  the  boundary  condition 
f(R,R)  =  gQ (R,0)  . 


(5.17) 
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,  £ 

This  shows  that  to  compute  E[z(0)|Y‘]  we  need  only  to  solve  simultaneously 
(5.16)  ar.d  the  vibrating  membrane  equation  (2.27)  for  g  (*,*).  After  dividing 

v. 

the  interval  [0,3]  into  N  subintervals  and  discretizing  these  equations,  this 
procedure  requires  only  0(N-)  operations  whereas  the  method  proposed  in  [9],  [10] 
would  require  C(N^)  operations.  In  some  sense,  the  relations  (5.16)  and  (2.27) 
constitute  a  generalization  of  the  Levinson  recursions  [23],  [24]:  g  (•,•)  and 
f{*,»)  correspond  respectively  to  the  forwards  and  backwards  predictors  of  [24], 
The  difference  is  that  here  the  domain  of  observation  grows  radially. 

The  result  described  above  can  be  used  as  follows:  assume  that  some 
observations  of  an  isotropic  random  field  are  given  over  a  large  domain,  e.g. 
the  whole  plane.  In  this  case,  we  could  construct  an  optimum  smoothing  filter 
based  on  all  observations,  but  to  estimate  the  field  everywhere,  a  very  large 
amount  of  computations  would  be  required.  An  approximate  smoothing  filter  can 
be  obtained  by  computing  f(R,*)  for  increasing  values  of  R  until  R  is  such  that 
the  mean- square  error 

g0(R,R)  =  E[(z(0)  -  E[z(0) 1yR])2] 

is  sufficiently  small.  Then  at  every  point  (r,9)  we  can  construct  an  estimate 
of  z(r,9)  by  using  the  filter  f(R,*)  and  the  observations  located  in  a  radius 
R  of  (r,9).  In  most  cases,  R  will  be  much  smaller  than  the  size  of  the  domain 
of  observation,  but  the  estimate  z  (r,9)  obtained  by  the  previous  procedure  will 
not  be  very  different  from  the  optimum  estimate  based  on  all  observations. 

VI.  Conclusion 

In  this  paper  we  have  obtained  a  set  of  vibrating  membrane  equations  for 
estimating  a  two-dimensional  isotropic  random  field  given  some  observations  of 
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this  field  over  a  finite  disk.  These  equations  depend  recursively  on  the  radius 
of  the  observation  disk  and  they  provide  an  efficient  method  for  constructing  the 
random  field  estimates.  They  also  generalize  the  vibrating  string  equations  of 
[1]  -  [5]  and  the  Levinson  recursions  [23]  which  had  been  introduced  to  solve  the 
corresponding  problem  in  one  dimension.  By  mapping  this  estimation  problem  to 
the  spectral  domain,  we  have  also  shown  that  the  estimation  procedure  described 
here  is  related  to  the  Gelfand-Levitan  solution  of  the  inverse  scattering 
problem  of  quantum  mechanics  for  circularly  symmetric  potentials. 

The  analysis  developed  in  this  paper  can  be  extended  easily  to  higher 
dimensional  isotropic  rancor,  fields  provided  that  instead  of  expanding  the 
random  field  in  Fourier  series,  we  expand  it  in  spherical  harmonics  as  shown 
in  [16] .  The  isotropy  of  the  random  field  has  been  reflected  here  by  the  fact 
that  the  vibrating  membrane  equations  are  circularly  symmetric.  If  the 
observed  random  field  is  only  homogeneous,  it  is  not  clear  whether  some  r.on- 
circularly  symmetric  membrane  equations  could  be  derived  to  solve  the  linear 
estimation  problem.  Mote  however  that  this  type  of  equation  appears  in  inverse 
scattering  studies  [14]  for  noncircular ly  symmetric  potentials.  Finally,  we  note 
that  the  connection  between  linear  estimation  and  inverse  scattering  problems 
has  only  been  discussed  briefly.  A  more  detailed  study  of  this  relation  will 
be  given  in  [17]  (see  also  [6]  -  [8]).  It  seems  likely  that  the  study  of  the 
connections  between  estimation  theory  and  inverse  scattering  problems  would  be 
beneficial  to  both  fields:  note  for  example  that  the  Kalman  filter  could  be  a 
valuable  tool  to  solve  problems  such  as  those  discussed  in  [25] . 

Appendix  A 

Proof  of  Theorem  1 

To  show  that  gn(*,  *)  is  twice  differentiable,  we  note  that  since  k(0  is 
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twice  differentiable,  i  (*,*)  is  also  twice  differentiable.  Then  (2.26)  can 

n 

be  rewritten  as 


-=-  a  (R,r i  =  :  (R,r)  -  •  g  (R,s)  ?  (s,r)  sds  (A.i) 

2"  'n  n  J  ^n  n 

where  the  right-hand  side  of  (A.I)  is  once,  then  twice  differentiable,  so  that 

g  (*,•)  is  twice  differentiable. 

Jn 

Now,  apply  the  displacement  operator  6  to  (2.26).  By  using  the  dis- 

_n 

placement  property  (2.13)  of  d  ,  we  obtain 


0  =  |  ( (A  (R)  q  (R, s) )  6  (s,r)  -  g  (R,s)  (1  (r)  $  (s,r))  sds 

_  -n  'n  n  n  -n  n 


(—  (Rg  (R, R) )  +^4"  <H, s) )  i  )  i>  <R,r) 

dR  n  n  n 


3R  n 


(A. 2) 


Rg_n(R,R)  ?n  (R,r)  +  ~  &n  9n<R,r)  . 


3v  taking  again  into  account  the  displacement  property  of  in»  and  integrating 
by  parts,  one  gets 


g  (R, s)  (A  (r)  (s,r))  sds 

n  -n  n 


I  (A  (r)  g  (r , s) )  $  (s,r)  sds  +  sg  (R,s)  S-  4>  (s ,r)  S  R 
J0  n  n  n  n  ds  n  !  s=0 


-  s  5—  g„(R,s)  <J>  (s,r) 
o  s  n  n 


so  that  (A. 2)  becomes 


0  =  ;  (<5  g  (R,  s) )  4>  (s,r)  sds  +  —  5  g  (R,r) 

1  —  t-*  n  n  jTjj  —  n  *1 


(R)  <h  (R,r) 


VRJ  ’n 
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where  a  (R)  is  given  bv  (2.28).  This  shows  that  5  g  (R,r)  is  the  solution 
*n  -n  n 

of  the  Fredholm  equation  (A. 3).  This  equation  is  identical  to  the  one  that 

would  be  obtained  bv  multiplying  (2.26)  by  a  (R) .  Therefore,  by  linearitv, 

~n 

and  by  noting  that  the  solution  of  (2.26)  is  unique,  we  have 


<5  g  (R,r)  =  q  (R)  g  (R,r) 
-n  =n  n  3n 


(A. 4) 


To  obtain  the  boundary  conditions  (2.29),  use  (2.26)  and  observe  the  $(•,•) 
satisfies  some  similar  conditions. 

Proof  of  Theorem  4 


By  operating  respectively  with  -  y-  and  on  the  integral 


3S 


equations  satisfied  by  g^  and  g  we  obtain 


Vr's)  =  j0(4 ' r}  9n(r'u)  VU's)  udu 


(A. 5) 


and 


+  ^  ^n(r,s)  +  rgn(r'r)  Vr's) 


, 9  (n+1) .  .  ,  .  fr  ,  .  , 9  (n+1) .  x  ,  . 

(__  +  - )  (j)  (r,s)  =  i  g  ,  (r,u)  (r—  +  - )  <t>  (u,s)  udu 

9s  s  n+1  J  n+1  9s  s  n+1 


1,9.  (n+1) ,  .  ,  x 

+  r—  (m—  +  - )  g  (r,s) 

2m  9  s  s  n+1 


(A. 6) 


Then,  if  we  take  into  account  the  property  (4.1a)  of  and  ^n+1  and  integrate 
by  parts,  we  get 


f r  ,  .  ,  9  (n+1) .  *  , 

j  gn+1<r,u)  (^  +  —  >  <Wu's)  udu 


(A. 7) 


((4+  ■^Ti)  gn+l(r'u))  VU'S>  udu  "  rgn+l(r'r)  Vr's)  * 


Now,  denote 


bn  -?  ( n+ 1  ) 

c  (r,s)  =  -  7)  g  (r , s)  + 

n  a  r  r  n  c?s  s 


Wr's) 


(A. 8) 


When  (A. 7)  is  substituted  inside  (A. 6) ,  by  adding  the  resulting  equation 

to  (A. 5)  we  find  that  c  (*,*)  satisfies  the  intearal  actuation 

n 

l*r 

0=1  c  (r,u)  <j>  (u,s)  udu  +  —  c  (r,s)  +  a  (r)  i>  ( r,s )  (A. 9) 

J  o  n  n  2"  n  nn 

where  a_^(r)  is  given  by  (4.4).  This  equation  can  be  viewed  as  obtained 
from  (2.26)  by  multiplication  by  -a^tr).  Therefore,  by  linearity,  and  since 
the  solution  of  (2.26)  is  unique,  we  have 


c  (r,  s)  =  -  a  (r)  g  (r,s) 
n  n  n 

The  identity  (4.3b)  can  be  derived  similarly. 


(A. 10) 
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